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We study the phase diagram of the superconducting vortex system in layered high-temperature 
superconductors in the presence of a magnetic field perpendicular to the layers and of random atomic 
scale point pinning centers. We consider the highly anisotropic limit where the pancake vortices on 
different layer are coupled only by their electromagnetic interaction. The free energy of the vortex 
system is then represented as a Ramakrishnan-Yussouff free energy functional of the time averaged 
vortex density. We numerically minimize this functional and examine the properties of the resulting 
phases. We find that, in the temperature (T) - pinning strength (s) plane at constant magnetic 
induction, the equilibrium phase at low T and s is a Bragg glass. As one increases s or T a first 
order phase transition occurs to another phase that we characterize as a pinned vortex liquid. The 
weakly pinned vortex liquid obtained for high T and small s smoothly crosses over to the strongly 
pinned vortex liquid as T is decreased or s increased - we do not find evidence for the existence, 
in thermodynamic equilibrium, of a distinct vortex glass phase in the range of pinning parameters 
considered here. We present results for the density correlation functions, the density and defect 
distributions, and the local field distribution accessible via /iSR experiments. These results are 
compared with those of existing theoretical, numerical and experimental studies. 

PACS numbers: 74.25.Qt, 74.72.Hs, 74.25.Ha, 74.78.Bz 



I. INTRODUCTION 

Equilibrium and transport properties of the mixed 
phase of highly anisotropic, layered, high-temperature 
superconductors (HTSCs) in the presence of various 
kinds of pinning have been studied extensively for over 
one decadeji Random pinning destroys the long-range 
translational order of the Abrikosov lattice and leads to 
a variety of low-temperature glassy phases. In systems 
with random point pinning, the existence of a topologi- 
cally ordered Bragg glass (BrG) phase with quasi-long- 
range translational order at low temperatures, low fields 
and weak pinning is now well-established both theoret- 
icallySi^ and experimentally^. The BrG is expected to 
melt into a vortex liquid (VL) via a first-order transition 
as the temperature is increased. The possibility of oc- 
currence of an amorphous vortex glass (VG) phase (with 
nonlinear voltage-current characteristics and vanishing 
resistance in the zero-current limit) at low temperatures 
in systems with strong pinning (or at high magnetic fields 
where the effects of pinning disorder are enhanced) was 
suggested several years age^-^. However, in spite of ex- 
tensive investigation over almost twenty years, whether a 
true VG phase (that is, an amorphous glassy phase ther- 
modynamically distinct from the high-temperature VL) 
exists in systems with uncorrelated point pinning remains 
controversial. 

Evidence for the theoretically expected scaling behav- 
ior of the current-voltage characteristics near a contin- 
uous VL-VG transition was obtained in several early 



experiments^ The validity of the conclusions drawn from 
these experiments has been questioned in latter studies^ 
although a recent experimental study^ has presented in- 
triguing thermodynamic evidence for the existence of a 
VG phase. Numerical studiesiS of a "gauge-glass" model 
believed to be in the same universality class as vortex 
systems with random point pinning suggest that a finite- 
temperature VG phase does not exist in three dimen- 
sions if electromagnetic screening of the interaction be- 
tween vortices is taken into account. Other numerical 
studiesiiiiSii^ of models in which the interaction between 
vortices is not screened indicate the presence of a thermo- 
dynamically stable VG phase at low temperatures if the 
pinning is strong. One of these^^ also presents evidence 
for the existence of a "vortex slush" phase with short- 
range translational order (i.e. a translational correlation 
length significantly longer than that in the VL phase) but 
no long-range phase coherence. However, the existence 
of this phase has been questioned in later workii4 A vari- 
ety of interesting "glassy" behavior has been experimen- 
tally observediS^ near the first-order melting transition 
of the BrG phase of superconductors (both conventional 
and HTSCs) with point pinning. It is possible tha t^^i^^ 
these observations can be understood by assuming that 
the melting of the BrG phase occurs in two steps: the 
BrG first transforms into a "multidomain" glassy phase 
that then melts into the usual VL at a slightly higher 
temperature. 

It is clear from this brief summary of the many ex- 
isting theoretical and experimental results that several 
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important issues about the structure and phase behavior 
of vortex matter in the presence of random point pinning 
are still unresolved. 

In this paper we present the results of a numerical 
study of the structural and thermodynamic properties 
of the system of vortices in a highly anisotropic layered 
superconductor in the presence of random point pinning. 
The magnetic field is perpendicular to the superconduct- 
ing layers and the vortex lines can then be considered as 
stacks of pancake vortices residing on the layers. These 
pancake vortices constitute a system of point-like objects 
that interact among themselves and also with an "ex- 
ternal" potential arising from a small concentration of 
randomly placed pinning centers (atomic scale defects). 
The positions of these defects on different layers are as- 
sumed to be completely uncorrelated. In this highly 
anisotropic case, the Josephson coupling between layers 
is vanishingly small, and pancake vortices on different 
layers are coupled only through the electromagnetic in- 
teraction. Under these circumstances, one can use (as 
in our recent studies^^'^^ of vortex systems in the pres- 
ence of columnar pinning) a spatially discretized version 
of the Ramakrishnan-Yussouff free energy functional^^ 
for a system of pancake vortices. In this formalism, the 
free energy of the system is expressed as a functional of 
the time-averaged local density of vortices. The differ- 
ent local minima of this functional can then be found 
by starting the numerical minimization process with ap- 
propriate initial conditions. These local minima, whose 
density configurations are obtained in our calculations, 
can then be identified with different phases of the system 
in our mean field description: the nature of the phase cor- 
responding to a particular local minimum is determined 
from a detailed analysis of the density distribution and 
the correlation functions at that minimum. A first-order 
transition between two phases, represented by two differ- 
ent locally stable minima of the free energy functional, is 
signaled by a crossing of the free energies of these min- 
ima as some control parameter (e.g. the temperature) is 
varied. 

Using this method, we have analyzed the phase dia- 
gram of the system as a function of temperature and 
pinning strength for a fixed value of the magnetic in- 
duction, which determines the areal density of pancake 
vortices. For low temperatures and relatively small val- 
ues of the pinning strength, we find a nearly crystalline 
minimum that exhibits all the properties expected for a 
BrG phase (the "BrG minimum"). This BrG phase is 
the thermodynamically stable one (has the lowest free 
energy) if the temperature is sufficiently small and the 
pinning sufficiently weak. It becomes unstable very soon 
after either of these parameters is increased beyond the 
region where this is the equilibrium phase. We find an- 
other local minimum at which the density distribution 
in each layer is strongly disordered (amorphous) and the 
densities in different layers are weakly correlated. This 
minimum is at least locally stable for the range of tem- 
perature and pinning strength we have considered. If the 



temperature is high and the pinning strength low, then 
the density distribution at this minimum is weakly in- 
homogeneous, showing the characteristics expected for a 
weakly pinned VL. As the temperature is reduced, or the 
pinning strength is increased, this minimum evolves con- 
tinuously into a state with strong inhomogeneity in the 
distribution of the time-averaged local density. This state 
represents a system of pancake vortices strongly localized 
at random positions because of the random pinning po- 
tential. This state exhibits characteristics of a glass, but 
we prefer to call it a strongly pinned VL because it is 
continuously connected to the weakly pinned VL found 
for high temperatures and weak pinning - we do not find 
any indication of a phase transition between the liquid- 
like and glass-like states. 

As the temperature is increased at a fixed, relatively 
small value of the pinning strength, or the pinning 
strength is increased at a sufficiently low temperature, 
the free energy of the BrG minimum crosses that of 
the VL minimum, signaling a first-order phase transition 
from the BrG to the pinned VL phase. For very weak 
pinning, the temperature at which this transition occurs 
is nearly independent of the pinning strength: it is very 
close to the melting temperature of the unpinned vortex 
lattice. As the pinning strength is increased, the transi- 
tion temperature decreases and eventually, the transition 
line in the temperature - pinning strength plane exhibits 
a trend towards becoming parallel to the temperature 
axis. Thus, the BrG phase is thermodynamically stable 
only if the pinning strength is lower than a critical value. 
The shape of the boundary of the BrG phase is similar to 
that found in existing numerical simulation studies^^*^ 
As mentioned above, the local minimum corresponding 
to the VL phase remains locally stable as the transi- 
tion line to the BrG phase is crossed either by reducing 
the temperature or by decreasing the pinning strength. 
On the other hand, the BrG minimum becomes unstable 
(and the minimization converges to a minimum of the 
VL type) as the temperature or the pinning strength is 
increased slightly beyond the point at which the tran- 
sition to the VL phase occurs. This behavior is consis- 
tent with experimental results2i that indicate that the 
superheated BrG phase in NbSe2 exhibits a spinodal in- 
stability, whereas there is no limit of metastability of the 
supercooled disordered state. We have also carried out a 
careful search for the presence of local minima of the free 
energy that may correspond to the "vortex slush" phasei^ 
or the polycrystalline "domain glass" phasei^ predicted 
in some earlier studies. As we will explain below, we did 
not find any equilibrium free energy minimum that could 
be interpreted as either one of these phases although we 
have found some indications that they may occur under 
nonequilibrium conditions. 

We have used our results for the vortex density distri- 
bution at the free energy minima to calculate the distri- 
bution of the local (microscopic) magnetic induction in 
the sample. This is important because the distribution 
of the local magnetic induction is measured in muon spin 
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rotation (/^SR) experiments, and phase transitions in the 
vortex system give rise to characteristic changes in the 
width and shape of this distribution. Thus, /iSR experi- 
ments have been widely used^^ to study phase transitions 
in the mixed state of HTSCs and conventional type-II su- 
perconductors. Our results show a relatively broad and 
asymmetric distribution of the microscopic magnetic in- 
duction in the BrG phase, and a narrow and much more 
symmetric distribution in the VL phase. These features 
are in good agreement with available /iSR data. 

The rest of this paper is organized as follows. In sec- 
tion ^ we define the model we consider and describe 
the numerical method used in our study. The results 
of our calculations are described in detail and compared 
with existing experimental and numerical results in sec- 
tion mil Section IIVI contains a summary of our main 
results and a few concluding remarks. 

II. MODEL AND METHODS 

We consider a layered superconductor in a magnetic 
field perpendicular to the layers. We assume that the 
material is strongly anisotropic so that the Josephson 
coupling between layers is vanishingly small. The pan- 
cake vortices on different layers are then coupled only 



through their electromagnetic interaction. In this limit, 
vortices in the same layer repel each other with a po- 
tential logarithmic in their separation, while there is a 
much weaker attractive interaction between vortices on 
different layers which falls off exponentially with layer 
separation and varies logarithmically with the separa- 
tion in the layer plane. This model for the vortex sys- 
tem was used in our earlier studiesiiiS of the effects 
of columnar pinning. As discussed in detail in Ref. ITtI 
many theoretical and experimental studies indicate that 
this model is appropriate for describing the properties 
of vortex matter in highly anisotropic materials such 
as BSCCO (Bi2Sr2CaCu208+x) in a large region of the 
field-temperature plane, although this conclusion has 
been questioned in a recent theoretical analysis^i. 

The free energy functional of such a system of pancake 
vortices lying on the superconducting layers and inter- 
acting with a time-independent pinning potential can be 
written in the form: 

F[p]=Fry[p\+FM (2.1) 

The first term in the right side of Eq. H2.1|l is the free 
energy of the vortex system in the absence of pinning, 
while the second includes the pinning effects. The first 
term is of the Ramakrishnan-Yussouff form^S: 



f3FRY[p]^Yl J rfr{p„(r)ln(p„(r)/po)-<5pn(r)}-(l/2)^^ J dr J dr'C^„(|r - r'|)^p„(r)(5p„(r'), (2.2) 



where (3 is the inverse temperature. We have defined 
Spnijc) = Pni^) — Po as the deviation of Pnii^), the time- 
averaged areal density of pancake vortices at point r 
on layer n, from pQ, the density of the uniform liquid 
{po = B/^o where B is the magnetic induction and $o 
is the flux quantum) and we have taken our zero of the 
free energy at its uniform liquid value. The integrations 
are two-dimensional and the sums are over all layers. The 
function Cm„(r), which depends on the layer separation 
|m — n\ and the separation r in the layer plane, is the 
direct pair correlation function of the layered vortex liq- 
uid'^^ at density pQ. This static correlation function con- 
tains all the required information about the interactions 
in the system. We have used here the Cmn{r) obtained 
from a calculation^^ via the hypernetted chain approxi- 
mation^ for parameter values appropriate for BSCCO, 
as in Ref. d = ISA for the interlayer distance, and 
a two-fluid form for the temperature dependence of the 
penetration depth A(r) with A(0) = 1500A. 



The second term in Ea. (|2.1|) represents the effects of 



pinning and can be written as 



^pW=E / rfrKf(r)[Pn(r) 



Pol 



(2.3) 



where V^{r) is the pinning potential at point r on layer 
n. The pinning potential is assumed to be produced by 
random atomic scale point defects. The positions of the 
defects on different layers are assumed to be completely 
uncorrelated. Each defect is assumed to produce a pin- 
ning potential of depth vq and having a short range which 
we will take (see below) as of the order of the spacing h 
of the computational mesh in our numerical calculation. 
A randomly chosen small fraction, c, of the unit cells of 
the underlying crystal lattice, taken to be a square lat- 
tice with spacing dg = AA, are assumed to be occupied 
by pinning defects. 

If the time averaged vortex density is the same on 
all layers (as would be the case for columnar pins per- 
pendicular to the layers), the free energy could be writ- 
ten in an effectively two-dimensional formiii^ with the 
quantity C{r) = pl^ying the role of the 

two-dimensional direct correlation function. But in the 
present case, with a pinning potential that varies from 
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one layer to another, the full three dimensional problem 
has to be considered. The problem is therefore compu- 
tationally much more demanding. To numerically min- 
imize the free energy of the system, we discretize the 
density variable by introducing a layered triangular com- 
putational lattice containing N'^ sites on each layer. The 
total number of computational mesh points is thus N'^N^ 
where is the number of layers considered. We use pe- 
riodic boundary conditions in all three directions. On 
the sites of this lattice we define discretized density vari- 
ables Pn^i = Pn(ji)AQ, where Pn(ri) is the density at mesh 
point i on layer n and Aq the area of the unit cell of the 
in-plane computational lattice. A numerical procedure 
developed earlieipii^ is used to find local minima of the 
free energy written as a function of the variables {pn,i}- 

We use the quantity Qq defined by ttpoUq = 1 as our 
unit of length. The spacing of the triangular computa- 
tional lattice is taken to be /i = a/16 where a = 1.998ao 
is the equilibrium lattice constant of the Abrikosov vor- 
tex lattice in the absence of pinning in the temperature 
range considered (see below) which is chiefiy determined 
by the melting temperature of the unpinned vortex lat- 
tice, namely ~ 18. 4K^^ at the field considered in this 
work {B = 0.2T). The total number of defects on each 
layer of our sample is given by Nd = V^{Nh)'^c/{2dl). 
In our discretized system, the values of the pinning po- 
tential Vn,i associated with mesh point i on the nth layer 
are determined in the following way: for each layer, we 
generate points distributed randomly and uniformly 
over the area of the layer. The potential 14,, i is then 
taken to be wo(™n.i — Nd/N'^) where mn,i is the number 
of such point in the ith computational cell on layer n, and 
Nd/N"^ is the average value of mn,i. We take c to be 0.01 
(1%) and measure vo in units of kg^ the Boltzmann con- 
stant, so that the parameter s = wo/fcs (with values given 
in Kelvins) measures the strength of the pinning poten- 
tial. With this parametrization, we have (yn,i) — 0, and 
j) = 2m{s/Tf5nn'5^j whcrc T is the temper- 
ature and the angular brackets denote the average over 
the random configuration of defects. 

To test the numerical method described above, we car- 
ried out a minimization of the free energy for a system in 
which a single strong pinning center was placed on one 
of the layers. The depth and the range of the poten- 
tial pro duced by the pinning center were chosen (as in 
Ref. [13, see Fig. 3 in that work) so as to ensure that 
it traps a single vortex at the temperature T = 20K 
(the unpinned vortex system is in the liquid state at this 
temperatureii). According to a well-known result in liq- 
uid state theorjiSi, the time-averaged local density p„(r) 
outside the range of the potential produced by the pin- 
ning center (here, n denotes the separation from the layer 
in which the pinning center is located, and r is the in- 
plane distance from the position of the pinning center) 
should be equal to Pog' n{f) where g' n(r) is the pair dis- 
tribution function of the unpinned vortex liquid, which 
can be obtained from the direct pair correlation function 
C„(r) used as input in our calculations. The results for 
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FIG. 1: (Color online). Plots of the local density p,i(r), 
normalized by the average density po, for n = ((red) circles) 
and n = 1 ((black) dots), as functions of r/ao. Here, n is the 
layer separation and r is the in-plane distance from a single 
strong pinning center that traps a pancake vortex. The results 
for the pair distribution function g'^if) are also shown for 
comparison ((red) dashed line: n = 0; (black) solid line: n — 
1). These results are for the vortex liquid state at temperature 
T = 20K and field B = 0.2T. 

Pn{f)l Po, obtained from a numerical minimization of the 
discretized free energy (with N = 128 and Nl — 32) 
are shown in Fig^ for n — Q and n — 1 (for n — Q, 
the results for pn [r) / po are shown for values of r outside 
the range of the potential produced by the pinning cen- 
ter at r = 0). The plots also show the values of g'„{r) 
obtained from the C„(r) used as input. The agreement 
between the two sets of results is quite good (the small 
differences may be attributed to the fact that the peak 
of the local density at the pinning center at r = is not 
a (5- function) . These results confirm that the numerical 
method used in our calculations correctly reproduces the 
strong intra-layer correlations, as well as the much weaker 
inter-layer correlations of local density of the system of 
pancake vortices. 

III. RESULTS 

We now present our results. The values of the temper- 
ature T and pin strength parameter s are indicated in 
each case. The field is always B = 0.2T. The number of 
layers is Nl = 128 and the transverse computational lat- 
tice size N = 256. With the value of the computational 
lattice constant h being h = a/16 where a = 1.988ao ifM^ 
the equilibrium lattice constant of the vortex lattice, the 
number of pancake vortices in each layer is = 256. 
The absence of finite size effects has been verified by 
checking that these larger-size results are nearly identical 
to those obtained for smaller systems. Results have been 
obtained for five random pin configurations and averaged 
over this number (except as indicated) whenever appro- 
priate. This is a sufficient number as we find that the dif- 
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ferences between results for different pin configurations 
of the same strengtli and concentration arc negligible for 
the quantities of interest. 

We use several kinds of initial conditions to start the 
minimization process. In the first kind, we start with uni- 
form "liquid-like" conditions, with the variables pn^i all 
equal to their average value. The second kind is "crystal- 
like" : the initial p„^i variables are obtained, for a given 
random pin configuration, by minimizing the pinning en- 
ergy of the equilibrium crystalline configuration in the 
absence of pinning^^ for the value of B used here, with 
respect to the symmetry operations of the computational 
lattice. Thirdly, once a local minimum configuration has 
been obtained from the minimization procedure for cer- 
tain values of s and T, that configuration is used as the 
initial condition for a nearby value of either one of these 
variables. The first kind of initial conditions leads, af- 
ter minimization, to free energy minima which are, as we 
shall see below, disordered in a liquid-like or glass-like 
way. The second kind yields BrG states, provided how- 
ever that the values of T or s are not too high: if they 
are, the crystal-like BrG state is not stable and those ini- 
tial conditions lead eventually to a disordered state. The 
third kind leads to a state of the same general kind as 
the initial one, except when an instability boundary is 
crossed. 



A. Structure of free energy minima 

^From the set of values for the pn^i at each local min- 
imum, the nature of that minimum can be ascertained 
by evaluating the appropriate correlation functions and 
also by direct visualization of either the density variables 
or of the "vortex lattice" formed by the locations of the 
peaks of the local density. 

To visualize the vortex lattice, the local peak densities 
are extracted from the set of pn,i values by defining site 
i in layer n to be a peak location when it corresponds 
to a local maximum of the pn^i such that the value pn^i 
exceeds that of any pnj in the same plane and within a 
radius a/2 from i. The positions of such "vortex sites" 
in each layer can then be plotted. An example is given 
in Fig. 121 The results there are all at s = 10. OX and 
T = 17.0/^. In the first panel, a small (blue) dot denotes 
the position of a vortex site at any one of the 128 layers. 
The results plotted correspond to a state obtained start- 
ing with "crystal-like" initial conditions, as explained 
above. The vortex positions are clearly clustered into 
well-defined groups that form a triangular pattern. For 
such minima, we can define "vortex lines" by joining the 
vortex positions in each group on different layers. The 
(red) filled circles in the first panel indicate the average 
positions of these vortex lines. These average positions 
form a slightly distorted triangular lattice, with consider- 
able variation, however, from layer to layer. The typical 
size of the clusters of (blue) dots represents the degree 
of transverse layer-to-layer wandering of the vortex lines. 



30 r 

25 
20 
15 

o 

^ 10 

5 

30 
25 
20 

«° 15 

03 
5, 

10 
5 
0. 



* »■ *: f .» * «■ » * * * * * 1- -t 4' 

t §■**-* 
i -i $ t 4 t 4 f f V 

•- -f » S f s '* • f # 
• *»*»■*' »M f *#■!! -* * * * 
« -f ir t - f f .# 4 .$ Ih # -it I! K $ -f 
# f * * f " i *. jf #■ * f f .# # * » 
-« -4' .# » i| it « ^' i -i .4 '11 # 1^ 

f t 4t * t- * **-*■"* ir »# i t * 
.* * t t * S .* .* * 4 * * * # * * 

g ^ f 1^ # |i; 4^ ^ ^ ^ f 4 % 



10 



20 



x/a. 



30 



40 



50 



A4A^&&^A4A^44i(L44 
^^A^A444^^^44^44 
4&44^Z*444«4444AA 
^A44d4A4444fc4&^^ 

A£^44^444i«AA4-9AC\* 
4&«'^/to4^^^A4ji»4^z]i 
4AAA4^4^^^4A4*^^^ 





30 p 

25- 
20 

>. 

10- 
5 



10 



20 x/a„ 30 



40 



50 




5 .^ maV 

qIa ^ . . ^ 

10 20 , 30 

x/a^ 



40 



50 



FIG. 2: (Color online). Visualization of density structures 
at different kinds of free energy minima. The first two panels 
correspond to an ordered minimum. In tlie first, tlie (blue) 
dots indicate, for each of the 128 layers present, the positions 
of pancake vortices in that layer. The (red) filled circles de- 
note the average positions of "vortex lines" obtained by con- 
necting adjacent vortices on neighboring layers. In the second 
panel, the (red) full circles have the same meaning and the 
(black) triangles represent the vortex position in a randomly 
chosen layer. In the third panel, the (red) full circles and the 
(black) triangles represent the positions of pancake vortices 
on two adjacent layers in a disordered minimum. In all cases 
T = 17.0K and s = lO.OA'. 
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In the second panel, the same point is emphasized by dis- 
playing again the average positions as (red) filled circles, 
and the vortex positions at a randomly chosen layer as 
triangles: the degree of crystalline order in the average 
positions of the vortex lines is higher than that in the 
vortex positions on a typical layer - the random trans- 
verse displacements of the vortex positions on different 
layers from the ideal triangular lattice partly cancel out 
in the calculation of the average positions of the vortex 
lines. From these results and other evidence discussed 
below, we identify minima of this kind as representing 
the BrG phase. 

If one plots the vortex positions as in the first panel 
of Fig. El for a minimum obtained from liquid- like initial 
conditions at the same T and s, one obtains a nearly 
uniform space-filling plot: the positions of the pancake 
vortices on different layers are nearly uncorrelated at this 
minimum. This difference between the two kinds of local 
free energy minima is emphasized in the third panel of 
Fig. 121 where the vortex positions in two neighboring lay- 
ers are plotted for a disordered minimum obtained for the 
same pin configuration, T, and s as those for the plots 
in the first two panels. Disordered local minima of this 
kind are identified below as representing the VL phase. 

It is very useful to discuss also the difference between 
ordered and disordered states in terms of correlation 
functions. We do this in the next two figures. Consider 
first Fig. 13 There we plot the "two-dimensional" static 
structure factor S'(k) defined as: 



5(k) = |p(k, fc, =0)|V(7V,7Vi) 



(3.1) 



where p(k, fc^) is the discrete Fourier transform of the 
{Pn,i\ in terms of the wavevector (k, /cx) (the z-axis is 
normal to the layers). In the figure T = 17. OX and s = 
10. OX. In the first panel S'(k) is plotted for the ordered 
state (which is, as we shall see, the equilibrium one in 
this case) while in the other panel the same quantity is 
plotted for the disordered state obtained from uniform 
initial conditions, which is locally stable at these values 
of s and T. The difference between the two states can 
be readily seen: note the large difference between the 
vertical scales in the two plots and the Bragg-like peaks 
in the ordered case. 

Alternatively, one can look at real-space correlation 
functions obtained from the discrete Fourier transform 
of 5(k, fcz). We define g{r,n) as the angularly averaged 
correlation of the time-averaged local densities at two 
points separated by n layers and in-plane distance r. The 
same-plane correlation function is g{r) = g{r,n = 0). 
For n 7^ 0, g{r,n) represents correlations between the 
vortex densities in different planes. We normalize g{r,n) 
by Pq, so that it approaches unity at large values of r 
for all n. These correlation functions are different from 
the pair distribution function g' n(r) shown in Fig. ^ 
the g{r, n) represent correlations of time-averaged local 
densities, whereas (?'„(r) are the equal-time correlation 
functions of the instantaneous local density. In particu- 
lar, g(r^ n) equals exactly unity for all n and r for a ho- 
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FIG. 3: (color online). The static structure factor as defined 
in Eq. I|3.1^ plotted versus the two dimensional wavevector (in 
units of a^^) at T = IT.OK and s = 10. OX for the equilibrium 
ordered state under these conditions (first panel) and for the 
metastable disordered state (second panel). 



mogeneous vortex liquid in the absence of any pinning, 
while g' n{r) for a homogeneous liquid shows oscillations 
(seen in Fig. ^ that reflect the short-range correlations 
present in a liquid. In the presence of a pinning poten- 
tial that makes the time-averaged local density in the liq- 
uid phase inhomogeneous, the correlation function g{r, n) 
represents the spatial correlations of this pinning- induced 
inhomogeneity. 

Results for g{r,n) are presented in Fig. 21 In the first 
panel g{r), g{r, 1) and g{r, 10) at T = 17.0 and s = lO.OX 
for the liquid-like metastable state are shown. The de- 
cay of the correlations with n is clear: the vortex densi- 
ties on different layers are only weakly correlated in the 
liquid-like phase, in agreement with the results shown in 
the third panel of Fig. [3 The corresponding plot for 
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FIG. 4: (Color online) Real-space angularly averaged cor- 
relation functions (see text) plotted versus two-dimensional 
dimensionless in-plane distance at T = 17.0K. The first 
panel displays (top to bottom curves at small r) g{r,n) for 
n = (red), n = 1 (green), and n — 10 (blue) at s = lO.OK, 
in the disordered state. The second and third panels show 
g{r) = g{r,n = 0) for the disordered and ordered states, 
respectively, at s = 10. OA" and s — 12.0K, with the paler 
(green) curve corresponding to the higher s. 



the ordered state at the same values of T and s (not 
displayed) would show no decay with n except at very 
small values of n, indicating that the vortices on differ- 
ent layers are close to registry in this state. In particular, 
g{r = 0, n) is found to be essentially independent of n for 
3 < n < Nl/2 in the ordered phase, showing absolutely 
no sign of decay for large values of n. In the other two 
panels, the variation of g{r) with s is shown for both the 
disordered and the ordered states and two values of s. 
The two cases are again in stark contrast. For the disor- 
dered, liquid- like minimum (second panel), the amplitude 
of the damped oscillations of g{r) is much smaller than 
that for the ordered, crystal- like minimum (third panel). 
The amplitude of these oscillations increases with s for 
the liquid-like minimum because the pinning-induced in- 
homogeneity in the local density increases as the pinning 
gets stronger, whereas for g{r) in the ordered minimum, 
this amplitude decreases as s is increased, reflecting that 
stronger pinning leads to larger distortions of the crys- 
talline structure. One could similarly plot the variation 
of g{r) with temperature: in this case the peak heights 
decrease as T increases. All these features of g{r,n) for 
the liquid-like state are very similar to those found in an 
analytic study^ « of the effects of pinning disorder on the 
correlations of a layered vortex liquid. 

From the correlation functions, an appropriate order 
parameter can be extracted. A convenient definition is 
m = g{rm) — 1 where is the first nonzero value of r 
for which g{r) has a peak. An alternative definitioni^ 
in terms of 5'(k, kz = 0) yields very similar results. In 
the first panel of Fig. \S\ this quantity is plotted as a 
function of the temperature at constant s = lO.OK for 
both the ordered and (inset) the disordered state. It is, 
as it should be, much larger in the former case and it 
decreases with increasing T for both states. As we shall 
see below, the transition temperature at this value of s 
is near T — \7AAK. In the second panel of Fig. |S1 to is 
plotted versus s at constant T = 17.0K for the same two 
cases. The transition as s is increased at this T occurs 
at s = ll.lK. We see that to now increases with s in 
the disordered state while remaining much smaller than 
its values in the ordered state. In the ordered state, m 
decreases with s, as expected. 

It is also useful to visualize the defect structure of the 
phases by means of Voronoi plots of the vortex lattice. 
In these plots the number of neighbors of each vortex 
lattice site in a given layer is found via a Wigner-Seitz 
construction. One then plots with different symbols the 
lattice points having six neighbors (this would be all the 
sites for a defect-free lattice) and those having e.g. seven 
or five neighbors (which correspond to single disclina- 
tions: neighboring pairs of disclinations of opposite sign 
correspond to dislocations). The results are shown in 
Fig. El The first three panels are for the ordered state, 
at s = lO.OK, T = 17. OK (top left panel), s = lO.OK, 
T = 17 AK (top right panel), and s ^ 12.0K, T ^ 17 .QK 
(bottom left panel). The system is deep in the ordered 
phase for the first set of values of s, T, while the other two 
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FIG. 5: (Color online). First panel: The order parameter m, 
as defined in the text, plotted versus temperature at constant 
s = lO.OK for the ordered state (main plot) and the disor- 
dered state (inset). Second panel: The same quantity plotted 
versus pin strength s at constant T = 17. QK for the ordered 
state (main plot) and the disordered state (inset). The solid 
lines are guides to the eye. 



sets represent points close to the phase boundary (see be- 
low) between the ordered and disordered phases. In each 
case, the defect structure in the layer with the largest 
number of defects is plotted. One can see a few dislo- 
cations in each of the three plots, but these dislocations 
occur in tightly bound clusters with zero net Burgers vec- 
tor. The number of defects is small and increases only 
weakly with T for fixed s, and with s if T is fixed. One 
clearly has a rather well ordered crystal-like structure 
even for values of s and T close to the phase boundary. 
A feature to note in these plots is that all these structures 
are single-domain: domain boundaries which appear in 
Voronoi plots as lines of defect pairs are clearly ab- 
sent. By contrast, the last panel in this figure, which 
depicts the disordered state obtained for s = lO.QK and 
T — n.&K (the disordered state is the equilibrium one 
for these parameter values, see below), shows a very large 
number of randomly distributed defects. The structure 
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FIG. 6: (Color online) Voronoi plots of the vortex lattice. 
The (black) dots denote ordinary sites with six neighbors, the 
(red) circles sites with seven neighbors and the (blue) triangles 
sites with five neighbors. The first three panels are for the 
ordered state with s = lO.OiC, T = 17. QK (top left panel), 
s = 10. OA', T = 17. AK (top right panel), and s = 12.QK, 
T — 17. OK (bottom left panel). The last panel is for the 
disordered state at s = IQ.QK, T = 17. QK. The (green) filled 
circles and (magenta) squares in this plot represent sites with 
eight and four neighbors, respectively. 
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FIG. 7: (Color online). The bond-orientational correlation 
function gair) ((red) triangles) and the translational corre- 
lation function gair) ((black) circles) as defined in the text, 
plotted as functions of in-plane distance r for the ordered state 
at s = 5.0/C, T = n.QK (filled symbols), and s = lO.Oif, 
T = 17. OA" (empty symbols). The dashed lines are guides to 
the eye. 
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FIG. 8: (Color online). The average width D of the local 
density peaks in the ordered phase (as defined in the text), 
measured in units of ao- In the main plot, results are shown 
versus T at pin strength s — 10. K while in the inset D is 
plotted versus s at T — 17.0 



of this state may be classified as amorphous - there is no 
indication of the presence of grain boundaries separating 
crystalline regions with different orientations. 

The degree of orientational order in vortex configura- 
tions is most conveniently defined in terms of the bond- 
orientational correlation function ge(r) defined as: 

ge{r) = (^(r)^(O)) (3.2a) 

where, as before, r is a vector in the layer plane, the 
brackets (• • • ) denote an average over the choice of the 
origin in a particular layer, over all the layers in the com- 
putational sample, and over angles. The field ?/'(r) is: 

^(r) = — ^exp[6z0,(r)] (3.2b) 

where Oj{r) is the angle made by the bond connecting a 
vortex at r to its j-th in-plane neighbor and a fixed axis, 
and n„ is the number of in-plane neighbors of the vortex 
at r, as determined from a Voronoi construction. Sim- 
ilarly, one can define a "translational correlation func- 
tion" 5g('") of the vortex lattice by an equation identical 
to the right-hand side of Eq. H3.2a|l . but with the field ip 
replaced by 

V-gW =exp(iG-r). (3.3) 

Here G is one of the shortest nonzero two-dimensional 
reciprocal lattice vectors of the triangular vortex lattice 
in the absence of pinning. We average over the six equiv- 
alent G's and over all the layers. 

Examples of these correlation functions are plotted in 
Fig. [7| where results are shown for the ordered minima 
at T = 17. OK and two different values, 5.0K and lO.OK, 



of the pinning strength s. The function ge{r) is shown 
by the (red) triangles and gcir) by the (black) circles. 
In each case, the upper plot corresponds to the smaller 
value of s. Both these functions appear to saturate to 
constant values at large r, indicating the presence of long- 
range bond-orientational order and (nearly) long-range 
translational order (the length scales considered here are 
probably too short to show the power-law decay of gc (r) 
expected^'^ for large r in a BrG phase) for both values of 
s. From these results, and those shown in Fig.|21(top two 
panels). Fig. |21 (top panel). Fig. 0] (bottom panel), and 
Fig. El (top three panels), we conclude that the ordered 
minima may be identified as representing the BrG phase. 

The width of the local density peaks of the vortex lat- 
tice at free-energy minima quantifies the degree of ther- 
mal and disorder-induced wandering of the pancake vor- 
tices from their time-averaged positions. We have studied 
how the average width of the local density peaks depends 
on T and s. The width D„_i of a local density peak at 
computational mesh point i on layer n is calculated from 
the relation 



l^j Pn,j 

where r^j is the distance between mesh points i and j on 
the same layer, and the sum is over the mesh points 
on the same layer that lie inside a unit cell of the vortex 
lattice centered at mesh point i. A similar calculation 
was carried out in our earlier study^^ of the Lindemann 
parameter at the melting transition of the vortex lattice 
in the absence of pinning. Here, the widths of different 
local density peaks are not all the same, and we define 
the average width D as the average of Dn,i over all the 
local density peaks in the sample. The variation of D 
with T and s is shown in Fig. |S1 where we plot D in 
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units of ao versus T at constant s = lO.OK in the main 
plot and at constant T — 17.0K versus s in the inset. 
As expected, D increases with T for fixed s. The value 
of D at a fixed T increases with s - the local density 
peaks become more distorted and broader as the pinning 
strength is increased. Melting of the ordered state (see 
next subsection) occurs when D reaches a value of ap- 
proximately 0.51 (this is true in all cases, not just those 
shown in the figure). Using the value of the equilibrium 
lattice constant, a — 1.988ao, for B = Q.2T, we get the 
value of the Lindemann parameter, £ = D/a, to be ap- 
proximately 0.26. This is close to the value obtained for 
the system without pinning^-- , indicating that the Linde- 
mann parameter at melting is nearly independent of the 
disorder strength s. 

The root-mean-square (rms) lateral displacement of a 
pancake vortex from the average position of the vortex 
line to which it belongs contains another contribution, 
arising from the pinning-induced transverse wandering 
of the position of the local density peak that represents 
the average position of the pancake vortex. As men- 
tioned in connection with Fig. ^ in the ordered phase 
the local density peaks on different layers are approxi- 
mately in registry. However plots such as those in Fig. [3 
indicate that there is an appreciable variation in the po- 
sition of a vortex line as one moves from layer to layer 
- the typical width of the clusters of (blue) dots in the 
first panel of Fig. [3 provides a measure of this variation. 
We have calculated the T- and s-dependence of D' , the 
rms displacement of the average positions of the pancake 
vortices that form a vortex line from the average position 
of the vortex line, further averaged over all vortex lines 
in an ordered state. The T-dependence of this quan- 
tity is weak in the range of temperatures considered - 
it increases slightly as T is increased. As expected, D' 
increases as s is increased at fixed T: its value in units 
of ao at T = 17. OK is of course zero at s = 0, 0.16 
for s = 5.0K, 0.23 for s = 7.5K, 0.27 for s = lO.OK, 
and 0.31 for s = 12.0K. Similar results are obtained at 
other temperatures. Thus, D' increases almost linearly 
with s for relatively small values of s and the rate of in- 
crease becomes smaller as s is increased further. D' is 
typically substantially smaller than D for the values of 
s and T studied. If the layer-to-layer random displace- 
ments of the average positions of the pancake vortices 
from the average position of the vortex line to which 
they belong and the lateral displacements of the pan- 
cake vortices from their average positions are assumed to 
be statistically independent, then the rms displacement 
of individual vortices from the avenge position of the vor- 
tex line would be given by Dnet = [D"^ + -D'^)^/^. This 
would lead to values of Dnet that are ^ 15% higher than 
those of D for the largest value of s {— 12.0if ) considered 
here. As discussed in the next subsections, these results 
are important for assessing the validity of the Lindemann 
criterion for the melting transition of the BrG state. 

We end this subsection with a few comments on the 
structure of the disordered VL minima. As seen above, 



the arrangement of the vortices at these minima may 
be classified as amorphous, not polycrystalline. While 
early decoration experiments^ on layered high- Tc super- 
conductors showed evidence for an amorphous structure 
of the disordered phase, more recent experiments^ on 
NbSe2 have shown the occurrence of polycrystalline dis- 
ordered structures. Also, simulations'^^''^^ of the struc- 
ture of two-dimensional vortex systems in the presence of 
random point pinning have found polycrystalline struc- 
tures. Our recent worliiSiSSiS on layered superconduc- 
tors with columnar pins perpendicular to the layers also 
showed polycrystalline disordered states. It is therefore 
pertinent to inquire why polycrystalline structures are 
not found in the present calculation. While a complete 
answer to this question is not yet available, we believe 
that this is a consequence of the nature of the pinning 
potential considered here. In a sample with a small con- 
centration of columnar pins perpendicular to the layers, 
the pinning potential (which is the same in all layers) 
is large at a few points and zero in the other parts of 
the sample. This favors the formation of polycrystalline 
structures: crystallites can form in the large "interstitial" 
regions where the pinning potential is zero, and since such 
regions on different layers are in registry, there is no con- 
flict between the preferred arrangements of vortices on 
different layers. In contrast, the formation of local crys- 
talline regions is less likely in the situation considered 
here because the pinning potential is nonzero at almost 
every computational mesh point. Also, since the pinning 
potentials on different layers are now completely uncorre- 
lated, the formation of crystallites would be energetically 
favorable (if at all) in different regions on different layers, 
and such regions would, in general, not be in registry. So, 
the interlayer interaction that favors vortices on adjacent 
layers being in registry would act against the formation 
of crystalline patches. The last factor is, of course, ab- 
sent in two-dimensional systems and less important for 
surface layers probed in decoration experiments. This 
may explain the observation of polycrystalline structures 
in two-dimensional simulations and some decoration ex- 
periments. 



B. Phase diagram 

We have seen then that there are two kinds of local 
minima of the free energy in the {s,T) region we have 
studied. The ordered minimum is to be identified as a 
Bragg glass (BrG) while the disordered one is a vortex 
liquid (VL, whether it should be called vortex glass is 
discussed below). 

Where both kinds of minima are locally stable, the 
equilibrium state is found by comparing their free en- 
ergies, which is trivial in our procedure. The transi- 
tions then are found by looking for free energy cross- 
ings. Examples are shown in Fig. |3 where results are 
shown as functions of T and s. We find that in all cases 
the transition is of first order: as in Fig. |^ the deriva- 
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FIG. 9: (Color online). Examples of the free energy per vor- 
tex, in units of ksT, plotted versus T at constant s = lO.OJC' 
(main plot) and versus s at constant T — 17. OK (inset). The 
(red) plus signs are data for the BrG phase and the (green) 
crosses for the disordered phase. The data points are con- 
nected by straight lines. Statistical error resulting from pin 
configuration average would result in error bars smaller than 
the symbols. First-order transitions occur at the crossings. 

tive of the free energy at the transition is discontinu- 
ous. In general, the BrG phase becomes unstable slightly 
above the transition: if one warms up that phase (or in- 
creases s) beyond the melting point one converges, after 
a large number of iterations, to the disordered phase. 
The reverse is not true, it is possible to substantially 
supercool the VL phase: the VL minimum remains lo- 
cally stable in the whole (s,T) region considered here 
(0 < s < 12.5K, 16.0K <T < 21.0K). This is consistent 
with experimental resultsSi that suggest that the BrG 
phase in NbSe2 exhibits a spinodal instability as the tem- 
perature is increased above the transition point, whereas 
there is no limit of metastability of the high-temperature 
disordered state as it is cooled to temperatures lower than 
the transition temperature. This is different from the 
behavior found in our recent studjiiSi^SiS^ of the vortex 
system in the presence of a small concentration of ran- 
dom columnar pins - in that case, the high-temperature 
VL minimum becomes unstable as the temperature is 
decreased below the transition point. The absence of 
a spinodal instability of the VL minimum in the present 
case is probably related to the fact (discussed above) that 
the formation of large crystallites which would make the 
VL minimum unstable is unlikely. 

From free energy results such as those shown in Fig.|3 
for a variety of values of s and T, it is possible to map the 
phase diagram in the s — T plane. The results, averaged 
over all samples and at constant B ~ 0.2T, are shown in 
Fig. IIUI We can see there that, upon increasing s from 
zero, the transition line separating the VL from the BrG 
is nearly vertical, and that it is only beyond s « 5. OA' 
that the transition line bends, while still retaining a quite 
appreciable T dependence of the value of s at which the 
transition occurs. The transition line trends towards be- 
ing parallel to the T-axis at lower temperatures, indicat- 




4 Bragg Glass 



2 - 

L , , , 

16.5 17 17.5 18 18.5 

T(k) 

FIG. 10: (Color online). Phase diagram in the s — T plane at 
constant field B = 0.2T. The symbols denote data points for 
the free energy crossings where first-order phase transitions 
between BrG and VL phases are found. They are connected 
by solid segments which thus represent the locus of first-order 
transitions. 



ing that the BrG phase is thermodynamically stable only 
if the pinning strength is lower than a "critical" value 
which appears to be somewhat higher than 12.0K. Our 
result for the "critical" value of s is in good agreement 
with an analytic estimate^ of this quantity for layered 
superconductors in the highly anisotropic limit where the 
electromagnetic interaction between vortices on different 
layers dominates over the interaction due to Josephson 
coupling. The general behavior of the transition tem- 
perature as a function of s is in agreement with theo- 
retical predictionsSS*^. The shape of the phase bound- 
ary in Fig. 111)1 is also very similar to that found in other 
numerical studiesi^iSS of the BrG-VL transition. The 
well-known result^ that the effective strength of pinning 
disorder in a sample with a fixed level of microscopic dis- 
order increases with the magnetic field B suggests that 
the phase boundary in the B — T plane would have a 
shape similar to that in Fig. IIUI with the s-axis replaced 
by the _B-axis. This would be consistent with experi- 
mental results^. However, caution should be exercised in 
translating our phase diagram in the s — T plane to one in 
the B — T plane because the analogy between increasing 
s at fixed B and increasing B for fixed disorder is not 
exact. As mentioned above, we find that the value of D, 
the average width of local density peaks at vortex posi- 
tions as defined below Eq. H3.4() , remains nearly constant 
along the phase boundary. This supports approximate 
analytic calculations^MS, based on the assumption that 
the melting transition of the BrG phase occurs at a con- 
stant value of the Lindemann parameter C. However, 
our results also show that if the quantity Dnet, which 
measures the rms displacement of vortices from the av- 
erage position of the vortex line to which they belong 
(see discussion in the paragraph following that including 
Eqn. H3.4|l ). is used to define C, then it would not remain 
constant along the phase boundary - it would increase 
by about 15% as the disorder strength is increased from 
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zero to 12.0K. Since the vortex positions on different 
layers are close to registry in the BrG phase, while the 
interlayer correlations in the vortex positions are very 
weak in the VL phase (see Figs. |21 and the BrG-VL 
transition may be thought of as a "decoupling" transi- 
tion^ at which the correlation between vortex positions 
on different layers drops abruptly. 

We do not find evidence for a VG phase thermody- 
namically distinct from the high-temperature VL in the 
range of pin strengths and concentrations studied. The 
VL minimum obtained at relatively high temperatures 
evolves continuously as T decreases. The degree of local- 
ization of the vortices, measured by the typical heights of 
the local density peaks at the free-energy minimum, does 
increase as T is reduced, but the change is smooth. We 
found no evidence of a phase transition between a weakly 
pinned and a strongly pinned distinct VL states. This 
leads us to classify all disordered minima at low temper- 
atures and relatively large pinning strengths as strongly 
pinned VL, to distinguish them from the more liquid-like 
disordered minima obtained for higher T and smaller s. 
We have carried out various thermal "cycling" procedures 
(such as warming up the disordered minimum obtained 
by quenching a liquid-like initial configuration to a low 
temperature, and cooling the VL minimum obtained at 
a high temperature) for various values of s to look for lo- 
cal minima that are different from the BrG and VL ones. 
We did not find (see also below) any other minimum that 
could be classified as a VG that is structurally and ther- 
modynamically different from the VL. Thus, we conclude 
that our system does not exhibit a distinct equilibrium 
VG phase in the region studied. Since the interaction be- 
tween straight vortex lines is screened in the model con- 
sidered here, such a conclusion would be consistent with 
the expectation^- that a VG phase does nor exist in three 
dimensions if electromagnetic screening of inter-vortex 
interactions is taken into account. However, such a dras- 
tically general inference would be unwarranted. Also, 
our mean-field calculation that considers only the den- 
sity variables can not rule out the occurrence of a contin- 
uous VL-VG transition with subtle characteristics (such 
as the continuous transition suggested by the experimen- 
tal results of Ref. 9, for which structural signatures and 
the nature of the order parameter are particularly hard 
to identify), or one in which the phase of the supercon- 
ducting order parameter (not included in our treatment) 
plays a crucial role. 

We mentioned in Section that earlier studies have 
suggested the occurrence of other phases for the system 
considered here. These include a "vortex slush" phasoM 
predicted to exist at intermediate values of T just above 
the boundary of the the BrG phase in the s T plane, 
and a "domain glass" phasei^ expected between the BrG 
and VL phases for relatively small values of s. It has been 
shown-- that the vortex slush "phase" found in the sim- 
ulation of Ref. corresponds to vortex configurations 
nearly crystalline in each layer, but with the vortices 
on different layers not being in registry: the orientation 




FIG. 11: (Color online). Intermediate non-equilibrated 
multi-domain state obtained at s = 9.0K, T = 17.0K as 
explained in the text. As in Fig. |5| the average position of 
each vortex is represented by a dot. 



of the reciprocal lattice vectors that describe the crys- 
talline order on different layers varies from layer to layer. 
The domain glass phase^^ is supposed to consist of large 
crystalline domains. We have looked for the existence 
of free-energy minima that may represent such phases. 
In our attempts, we artificially made the VL minimum 
unstable towards the formation of crystalline structures 
by multiplying the C,„„(r) used as input in our calcu- 
lations by a factor of 1.05, which makes the structure 
factor SCkjkz) = 1/[1 — poCCk, k^)] negative for fc^ = 
and k close to the magnitude of the smallest reciprocal 
lattice vectors of the triangular vortex lattice. Our expec- 
tation was that the system would then converge to one 
of the other kinds of minima mentioned above, if they 
existed. In these calculations, we started with liquid-like 
initial conditions and quenched to within the region (low 
T, moderately large values of s) where these phases are 
expected. The minimization routine was run for a large 
number of iterations with the modified Cmn(?'), so that 
the instability of the VL minimum resulted, indeed, in 
the formation of locally crystalline structures on the lay- 
ers. Gonfigurations obtained at different stages of this 
run were then used as inputs in minimization runs in 
which we slowly eliminated the modification in Cmnir) 
so as to recover the original, physical interaction. In all 
such runs, the system converged to the BrG state but 
very slowly: the number of iterations it took for con- 
vergence was over twenty times larger than that in the 
ordinary case. The system lingered for a very large num- 
ber of iterations in a multi-domain state. The ultimate 
convergence was always to minima with the vortices on 
different layers in registry, even when the vortices were 
very much out of registry in the intermediate configu- 
ration. An example is shown in Fig. ^] at s = 9.0if , 
T — 17.0K where the vortex positions for an interme- 
diate state with multi-domain structure are shown in a 
plot similar to those in Fig.|21 The converged state, how- 
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FIG. 12: (Color online) The probability distribution for fields 
P{b) in terms of the dimensionless field 6 at T = ISAK and 
s = 5. OA", which is very near the melting point. The (red) 
solid lines are for the BrG state and the (green) blocks are for 
the VL state. 



ever, turned out to be a single domain BrG very similar 
to that shown in the first two panels of Fig. |21 

These results argue against the existence of the other 
phases suggested in earlier studies, at least in the region 
of pinning strength and concentration studied, as ther- 
modynamic equilibrium phases. However, that the sys- 
tem remains in a polycrystalline state for a large number 
of iterations in the minimization process suggests that 
such states may be stabilized through changes in the 
model parameters or may correspond to long-lived tran- 
sient situations. Further investigations of this possibility 
would be interesting. On the other hand, the reported 
observation of this phase in simulations might result from 
using a small system size or insufficient equilibration, or 
both. In our samples, the number of layers and the num- 
ber of vortices on each layer are substantially larger than 
those in the simulations of Ref. Also, the interlayer 
interaction in our model is rather weak, since we consider 
only the electromagnetic interaction between vortices on 
different layers. The system nevertheless converges to a 
BrG minimum with vortices on different layers in reg- 
istry, which implies that a vortex slush in which vortices 
form crystalline arrangements on the layers, but the ar- 
rangements on different layers are not in registry is rather 
unlikely to occur as a thermodynamic phase. 



C. Distribution of local magnetic induction 

Finally, we give a brief discussion of our results for the 
distribution of the local, microscopic, magnetic induction 
in the superconductor. This is a quantity of considerable 
experimental interest, as it can be measured directly in 
/iSR experiments^^. The probability distribution of the 
local magnetic induction is easily computed'^® from the 



vortex density distribution at the minima. We consider 
the the component of the magnetic induction perpen- 
dicular to the layers. The Fourier transform of this z- 
component due to a single pancake vortex at the origin 
is given byS 



l + A2p + A2/fc2' 



(3.5) 



where d is the spacing between the layers and A the pene- 
tration depth in the layer plane. If the typical time scales 
for the dynamics of density fluctuations in the state the 
vortex system is in are much shorter than the character- 
istic time scale of muon spin precession, the muons see 
a broadened (time-averaged) density distribution of the 
vortices. The time-averaged local magnetic field Bz is 
then obtained through a convolution of the local density 
at a free-energy minimum with the expression for bz given 
in Ea. (|3.5(l . If the dynamic time scales for such fluctu- 
ations are much longer than the precession time of the 
muon spins, then the convolution is to be done^S using 
a local density consisting of two-dimensional (5-functions 
located at the local density peaks at the free-energy min- 
imum. The time scales are beyond the scope of our time 
averaged calculation, but both methods give very similar 
results in any case. 

We have used these procedures to calculate the dis- 
tribution of Bz for different free-energy minima. In our 
discretized system, the convolution is done using discrete 
Fourier transforms, leading to the values of Bz at the 
computational mesh points. These values are binned to 
yield histograms for the distribution of Bz- The results 
are best represented in terms of the dimensionless field 
b = Bz/B. An example is given in Fig. El at a point on 
the melting curve. The two states clearly have different 
probability distributions P{b), with that corresponding 
to the VL state being narrower and much more symmet- 
ric. The distribution for the BrG state shows a long 
"tail" in the large-6 side. All these features are very 
similar to experimental results^^ ''^. The results shown 
in Fig. El were obtained using the locations of the vor- 
tices at the minima - the distributions obtained from the 
time-averaged local density have the same shape, but are 
somewhat narrower in the VL phase. 

The degree of asymmetry of P{b) can be conveniently 
characterized^^ by introducing a dimensionless parame- 
ter a defined as the ratio of the cube root of the third 
moment of P{b) to the square root of its second moment. 
This parameter changes discontinuously at melting. An 
example is shown in Fig. 1131 where a is plotted as a func- 
tion of T at constant s = lO.OK. Results are given for 
the equilibrium state at each temperature and averaged 
over four pin configurations. It turns out that for the VL 
state (but not for the BrG), the results for a vary appre- 
ciably from one pin configuration to another (the only 
quantity we have found that does so) and also depend- 
ing on which of the two methods (local density peaks or 
full density) is used. We have therefore averaged also 
over the two methods. The resulting error bars are still 
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FIG. 13: (Color online). The parameter a characterizing the 
degree of asymmetry of the probability distribution P{b) (see 
Figureinj plotted as a function of T at s = lO.OK. The result 
plotted at each temperature is the value for the equilibrium 
state at that temperature, averaged over pin configurations. 

larger in the VL state as the distribution in that case 
is close to being symmetric and rather small changes in 
the shape of P{b) lead, in that situation, to considerable 
changes in a. The results for a are quantitatively similar 
to those of /iSR experiments22, on the BrG-VL transition 
in BSCCO induced by increasing the magnetic field at a 
fixed temperature. A similar discontinuity in the value 
of a is also observed^ at the melting of the BrG phase 
upon increasing the temperature at fixed B. However, 
the value of a in the VL phase is found^ to be negative 
in that case. The reason for the negative value of a is not 
completely clear. It has been attributed to effects of sam- 
ple geometrjiSS and it has been argued^ that a negative 
a indicates the presence of a multi-domain state. Since 
effects of sample geometry are not present in our calcu- 
lations with periodic boundary conditions and we do not 
find any multi-domain state, it is perhaps not surprising 
that our calculations do not yield negative values of a. 



IV. CONCLUSIONS 

We have studied in this paper the phase diagram of the 
superconducting vortex lattice in a layered superconduc- 
tor, in the presence of both a magnetic field normal to 
the layers and of random atomic-scale point pinning cen- 
ters. We have considered the strongly anisotropic case 
where the interactions are of two-body form and the free 
energy can then be written as a functional of the time 
averaged vortex density. Numerical minimization of the 
discretized free energy then yields the pancake vortex 
density distribution at the minima. 

The values of the free energy at the minima can be 
plotted to locate the phase transitions. Our main result 
is the phase diagram (Fig. I10|l in the temperature ~ pin 



strength plane, where we find a line of first-order phase 
transitions between a Bragg glass and a disordered phase. 
The density distributions at the minima can then be an- 
alyzed to obtain any desired correlation functions. They 
can also be used to directly visualize the vortex lattice 
in real space and to obtain the distribution of the micro- 
scopic local magnetic induction in the superconductor. 
From an extensive analysis of such results we conclude 
that the disordered phase is to be identified as a pinned 
vortex liquid. While one can distinguish between strong 
and weak pinning regimes, the crossover between them 
is completely smooth. 

As discussed in detail in the previous section, our re- 
sults are consistent with those of experiments and simu- 
lations, the only possible exception being that we do not 
find, in the region of pinning strength and concentration 
studied, any sign of additional "vortex glass" , "vortex 
slush" or "domain glass" phases the existence of which 
has been made at least plausible by results obtained 
e.g. in Refs. Igl lTMT^ respectively. We have ourselves 
found^^''^^-^'^ a polycrystalline Bose glass phase in the case 
of columnar pinning through the use of the same meth- 
ods as in this work. Such phases may exist under other 
conditions, specifically for a much smaller concentration 
of considerably stronger defects which would favor poly- 
crystallinity much more than the situation studied here, 
where we have chosen a value of the concentration c of 
atomic defects which, although small in the atomic scale, 
is in a sense large with respect to the vortex lattice: there 
are many pinning sites per vortex lattice unit cell. The 
results for columnar pinning, where the concentration is 
much smaller but the pinning strength larger, make it 
possible to conjecture that there is a regime of larger s 
but smaller c where additional phases may exist. Further 
work in that direction is clearly worthwhile. 

In a paper**^ that appeared after this one was origi- 
nally submitted, it has been suggested, from an analysis 
of a disordered xy model in which electromagnetic in- 
teractions between vortices is neglected, that the BrG 
phase should become unstable (no matter how weak the 
pinning) to disorder in the direction normal to the lay- 
ers at a value of B smaller than that considered here. 
Our results for the correlation function g{0, n) as a func- 
tion of interlayer separation n in the ordered phase (dis- 
cussed in the text in connection with Fig.^, indicate no 
decay in this quantity as a function of n in the region 
3 < n < Nl/2. Further, as mentioned above, our phase 
diagram data quantitatively agree with analytic calcula- 
tions-^ in which the softening of the tilt modulus at large 
wavelengths was included. These results strongly argue 
for the actual stability of the BrG phase, in agreement 
with previous numerical studiesi^iS£ of the disordered xy 
model. However, even though our samples are already, 
as noted, much thicker than those considered in earlier 
numerical studie s"^'^° , the possibility that the BrG phase 
would disorder in the transverse direction over a length 
scale much larger than A can not be mathematically ruled 
out by our numerical results. 
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